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Abstract 

Programs that calculate observables in quantum chromo dynamics at next-to-leading order typ- 
ically suffer from the problem that, when considered as event generators, the events generated 
consist of partons rather than hadrons and just a few partons at that. Thus the events gener- 
ated are completely unrealistic. These programs would be much more useful if the few partons 
were turned into parton showers, in the style of standard Monte Carlo event generators. Then 
the parton showers could be given to one of the Monte Carlo event generators to produce hadron 
showers. We show how to generate parton showers related to the final state collinear singularities 
of the next-to-leading order perturbative calculation for the example of e + + e~ ^3 jets. The soft 
singularities are left for a separate publication. 
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I. INTRODUCTION 



Perturbation theory is a powerful tool for deriving the consequences of the standard model 
and its extensions in order to compare to experimental results. When strongly interacting 
particles are involved, perturbation theory becomes relevant when the experiment involves a 
short distance process, with a momentum scale Q that is large compared to 1 GeV. Then one 
can expand the theoretical quantities for the short distance part of the reaction in powers 
of a s (Q), which is small for large Q. Of course, the particles continue to interact at long 
distances, but for suitable types of experiments the long distance effects in the final state can 
be neglected while the long distance effects in the initial state can be factored into parton 
distribution functions that describe the distributions of partons in incoming hadrons. This 
works if the measurement of properties of the final state is infrared safe: the measurement 
is only weakly affected by long distance interactions in the final state. Examples include the 
inclusive production of very heavy particles and the production of (suitably defined) jets of 
light particles. 

The simplest sort of calculation for this purpose consists of calculating the hard process 
cross section at the lowest order in a s , call it af , at which it occurs. For example, for two-jet 
production in hadron-hadron collisions, B = 2 while for three-jet production in electron- 
positron annihilation B — 1. This kind of lowest order (LO) calculation is simple, but leaves 
out numerically significant contributions. Corrections from higher order graphs are often 
found to be 50% of the lowest order result. For this reason, one often performs a next- 
to- leading order (NLO) calculation, including terms proportional to af and ctf +1 . Then 
the estimated error from yet higher order terms (which are usually unknown) is typically 
smaller, say 10%. 

There is another kind of calculational technique available, that of the Monte Carlo event 
generator [1-3]. In this technique, a computer program generates "events," a list of particles 
(ir, p, p, . . . ) and their momenta Pi,P2,P3, ■ ■ ■ that constitute a final state that could arise 
from the given collision. Call the nth final state f n . The program is constructed in such a way 
that the probability of generating a final state / is approximately proportional to the physical 
probability of getting this same final state in the standard model or in whatever theory is 
to be tested. That is, the predicted value for an observable described by a measurement 
function S(f) is 

(S) = ^lts(fn)- (1) 

JV n=l 

Alternatively, the program may generate weights w n for the events, so that 

1 N 

(s) = -j:^ l s(f n ). (2) 

JV n=l 

The Monte Carlo event generators have proven to be of enormous value. Among their 
benefits is the following. Since they generate complete final states, one can use detector 
simulations to analyze the effects on the measurement of any departure of the detector from 
an ideal detector. Essentially, this amounts to replacing an idealized measurement function 
S(f) by something that is much closer to what an actual detector does. 

There are, of course, some weaknesses in Monte Carlo event generators, considered as 
tools to provide predictions of the standard model and its extensions. One is that the 
theory encoded in them is more than the standard model. They also contain models and 
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approximations for long and medium distance physics. Fortunately, the long and medium 
distance effects do not matter much in the case that the measurement function S(f) is 
infrared safe and involves only one large momentum scale Q. In this case, the place where 
the medium and long distance approximations are important is in assessing the effects of 
detector imperfections, which likely result in having an effective measurement function, 
<S e ff(/), including the detector effects, that departs somewhat from the ideal function S(f) 
and thus departs from being truly infrared safe. 

A more serious weakness of current Monte Carlo event generators is that they are based 
on leading order perturbation theory, so that a rather substantial theoretical uncertainty 
should be ascribed to the result if the result is considered as deduction from the standard 
model (or one of its extensions) with no other hypotheses included. 

Return now to the NLO calculations. They are typically computer programs that act 
as Monte Carlo generators of partonic events. They produce events / consisting of a few 
partons with specified momenta. For example, for three-jet production in electron-positron 
annihilation, there are three or four partons in the final state. Each event comes with a 
weight w n . Then the predicted value for an observable described by a measurement function 
S(f) is 

1 N 

<S> = T7E«'"SU)- ( 3 ) 

JV n=l 

Here S(f) is an idealized measurement function applied to the partonic final state. We 
see that the NLO programs are actually quite close in structure to the Monte Carlo event 
generators. 

One perceived weakness of the NLO Monte Carlo generators is that the weights are 
negative as well as positive. This feature arises because the programs do real quantum 
mechanics, in which an amplitude times a complex conjugate amplitude has a real part that 
can be either positive or negative. Negative weights do not, however, cause a problem when 
inserted into Eq. (3): computers can easily multiply negative numbers. 

The real weakness of most current NLO Monte Carlo generators is that they produce 
simulated final states that are not close to physical final states. What could be done to 
generate realistic final states? Following the example of the leading order event generators, 
which do generate realistic final states, one would use the final states in the NLO partonic 
generators to generate a parton shower from each outgoing parton. This gives a final state 
with lots of partons. Then one could use one of the highly successful algorithms of the 
leading order Monte Carlo event generators to turn the many partons into hadrons. There 
should not be a serious problem with the hadronization stage, since this is modeled as a long 
distance process that leaves infrared safe observables largely unchanged. The problem lies 
with the parton showers. Here a high energy parton splits into two daughter partons, which 
each split into two more partons. As this process continues, the virtualities of successive 
pairs of daughter partons gets smaller and smaller, representing splittings that happen at 
larger and larger distance scales. The late splittings leave infrared safe observables largely 
unchanged. However, the first splittings in a parton shower can involve large virtualities. 
They represent a mixture of long distance and short distance physics. Thus one must be 
careful that the showering does not reproduce some piece of short distance physics that was 
already included in the NLO calculation. 

Consider an example that illustrates these ideas, jet production in pp collisions with 
y^s = 14 TeV at the Large Hadron Collider now under construction. One can calculate the 
inclusive cross section da/dMjj to make two jets with a dijet mass Mjj. We may specify that 
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the jets are defined with the kr algorithm [4], that we integrate over the rapidities of the two 
jets in the range \yi\ < 1, 1 2/2 1 < 1, and that Mjj is in the range 0.5 TeV < Mjj < 5 TeV. 
There are computer programs currently available to perform this calculation at NLO [5- 
7]. There are no large logarithms in the theory in this kinematic region. Thus an NLO 
calculation should be reliable. If the experiment were to disagree with the theory beyond 
the experimental errors and the estimated errors from uncalculated NNLO contributions and 
from parton distributions, that would be a significant indication that some sort of physics 
not included in the standard model is at work. 

Although the theoretical picture just outlined seems bright, the available programs do 
not produce realistic final states. To illustrate this, imagine computing a cross section that is 
sensitive to the final state. For example, one might use a standard NLO program to calculate 
the cross section da/dM JJ dM 1 , where Mi is the mass of the jet with the larger rapidity. 
Our result would not be sensible in the region Mi <C Mjj. It would have a nonintegrable 
singularity at Mi = together with a term proportional to 5 (Mi) with an infinite negative 
coefficient. The integral over Mi of this highly singular function would give back da/dMjj, 
but the differential distribution would be highly nonphysical. In contrast, the standard 
Monte Carlo event generators can generate sensible results simultaneously for both da/dMjj 
and da/dMjj dMi. The only drawback is that the result for da/dMjj will be correct only 
to leading order in an expansion in powers of a s (Mjj). The problem addressed in this paper 
(using, however, an example from electron-positron annihilation) is how to modify such an 
NLO calculation so as to take advantage of the parton shower technology built into Monte 
Carlo event generators. We wish to get quantities like da/dMjj dM x approximately right 
while not losing the NLO accuracy of da/dMjj = J dMi da/dMjj dM t . 

Some technical language is helpful here in order to make these ideas more precise. Using 
the familiar concept of infrared safety (defined precisely in Sec. II), one can describe the 
cross section da/dMjj as an infrared safe two-jet cross section. Then the cross section 
da /dMjjdMi for Mi > is an infrared safe three-jet cross section: it is infrared safe and is 
nonzero only when the number of final state particles is three or more. With this language, 
what we wish to do is get two-jet cross sections right to NLO and at the same time get n-jet 
cross sections approximately right for n > 2. 

Recall that the shower Monte Carlo programs work by using an approximate version of the 
perturbative theory at an infinite number of orders in perturbation theory, approximating 
the behavior of the Feynman graphs near the singularities corresponding to collinear parton 
splitting or soft gluon emission. These singularities lead to large logarithms, two powers of 
log(Mi/Mjj) for each additional power of a s . The programs approximately add up a series 
of factors times logarithmic factors. After this, there is a phase that models how partons 
turn into hadrons. In this paper we omit the hadronization phase and use our own code 
for the parton showering. Our intention is that in the future only the stages of showering 
that directly connect to the hard scattering would be part of a program for NLO hard 
scattering with showers. The parts of the shower related to "softer" parton splittings, along 
with hadronization, would be performed by one of the standard Monte Carlo programs. We 
thus concentrate in this paper on the showering-hard-scattering interface and not on the 
showering itself. 

We can illustrate the same points for an NLO base calculation that has one more jet. 
An instructive example is the cross section discussed above, da j dMjj dM\. Here we do not 
integrate over Mi and we insist that M\jMjj not be small. This is technically an infrared 
safe three-jet cross section, as discussed earlier. In a three or four parton system that appears 
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in the calculation of this quantity at NLO, jet 1 consists of two or three partons and has mass 
Mi, while jet 2 consists of one or two partons. The NLO calculation (using [7]) should be 
reliable but would not have sensible final states. Thus if we let M 2 be the mass of whichever of 
the two measured jets has the smaller rapidity and tried to calculate da/dMjj dMi dM 2 , the 
result would not be even approximately correct for M 2 <C Mjj and M 2 <C Mi. One would 
then like to have a method that could get quantities like da/dMjj dMi dM 2 approximately 
right while not loosing the NLO accuracy of da/dMjjdMi. 

This discussion can serve to illustrate another point. Note that a calculation of 
da/dMjjdMi will have a problem if Mi/Mjj <C 1 because the result will contain large 
logarithms, namely logs of Mi/Mjj. The large logarithms are associated with the partons 
getting close to a two-jet configuration. One could then view the event as arising from a 
two-jet event in which one of the jets splits into two nearly collinear jets or one hard jet 
and one soft jet. A two-jet calculation with showers would get the effects of these large log- 
arithms approximately right. What one would then want is a method to merge the two-jet 
and three-jet calculations so that 1) a two-jet infrared safe observable is correctly calculated 
at order a^, i.e. NLO, 2) a three-jet infrared safe observable with suitable cuts so that the 
contributing events are well away from the two-jet region is correctly calculated at order a A s , 
i.e. NLO for three-jets, and 3) a three-jet infrared safe observable with cuts that allow con- 
tributions from close to the two-jet region is correctly calculated at order a A s with additional 
corrections proportional to a™ +4 times logs, with n > 1, that approximately account for the 
production of three jets by showering from a two-jet event. An investigation of this merging 
problem would certainly be desirable, but is beyond the scope of this paper. 

This paper concerns e + + e~ — > 3 jets. The goal is to see how to modify a NLO calculation 
by adding parton showers in such a way that when the generated events are used to calculate 
an infrared safe observable, the observable is calculated correctly at next-to-leading order. 
To do this, we build on our work in [8] , which sets up the NLO calculation in the Coulomb 
gauge. It is convenient to use a physical gauge for this purpose because it makes the analysis 
simpler. In this paper we explain the aspects of this problem that relate to the divergences 
that appear in perturbation theory when two massless partons become collinear. This part 
of the method is, we think, quite simple and easy to understand. The treatment of the 
divergences that appear when a gluon becomes soft is a bit more technical and will be given 
in a separate paper [9]. 

The algorithms that we propose apply to showers from final state massless partons. 
They are illustrated using the process e + + e~ ^3 jets in quantum chromodynamics and 
are embodied in computer code that performs calculations for this process [10]. We do not 
discuss showers from initial state massless partons. While our research program (starting 
with [8]) was underway, Frixione and Webber succeeded in matching showers to a NLO 
calculation for a problem with massless partons in the initial state (but no observed colored 
partons in the final state) [11]. Subsequently, Frixione, Webber, and Nason have extended 
this method to massive colored partons in the final state [12]. The general approach of 
[11, 12] is quite similar to that of the present paper with respect to collinear singularities. 
The principle difference is that we generate the first splittings of partons from an order 
af diagram using the functions given by the partonic self-energy diagrams in the Coulomb 
gauge, while [11, 12] uses the splitting functions of the Herwig Monte Carlo event generator. 

We discuss e + + — > 3 jets but not e + + e~ — > 2 jets. Of course, to calculate e + + — > 
2 jets at order a s with parton showering added would be quite a lot simpler than the present 
calculation for e + + e~ — > 3 jets at order a 2 s . However, merging the two calculations would 
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not be trivial and is not attempted here. We simply suppose that the present calculation is 
to be applied with a measurement function that gives zero for events that are close to the 
two-jet configuration (for example, for thrust near one). 

There has been other recent work on this problem [13]. That of Collins [14] is particularly 
instructive. This approach seeks both to put the parton showering calculation on a more 
rigorous basis and at the same time to match it to the hard scattering calculation, viewing 
the hard scattering as one end of a cascade of virtualities. Our goal is more limited. We 
take it as proven by experience that parton showers work well to make reasonably realistic 
final states. Our only concern is to make sure that when we add parton showers we do not 
lose the NLO accuracy of the result. 

II. NEXT-TO-LEADING ORDER CALCULATION 

We begin with a program that calculates an infrared safe observable in electron-positron 
annihilation at next-to-leading order. The program [8] uses a physical gauge, specifically 
the Coulomb gauge, and performs all of the integrals, including the virtual loop integrals, 
numerically [15-17]. It is not really necessary to perform the virtual loop integrals numeri- 
cally. It is a lot easier, but if one wanted to perform these integrals by hand ahead of time 
and insert the analytical results into the program, the calculation could still proceed. It 
is also not necessary to use a physical gauge. However, the use of a physical gauge makes 
the method conceptually simpler. In particular, collinear singularities are isolated in very 
simple Feynman subdiagrams, the cut and uncut two-point functions. If one were to use the 
Feynman gauge, the collinear singularities would be spread over many kinds of graphs and 
one would need another layer of analysis to organize them into a simple form. The Coulomb 
gauge, in particular, is well suited for calculations of electron-positron annihilation, in which 
the physics is approximately rotationally invariant in the electron-positron cm. frame. The 
time-like axial gauge has this same good property, but is not well suited for our algorithm for 
a technical reason: it introduces singularities in the complex energy plane for loop integrals, 
which would make the evaluation of the energy integrals in our algorithm more complicated. 

In order to present a method for adding showers to the NLO program, we need to review 
a little about how the NLO calculation works. We present the barest outline. Details may 
be found in [8, 15-17]. In the NLO program, the observable is expressed in the form 

J = H -7 / dpi ■ ■ ■ dp n 5 (^2 Pi) k(y /*X 

xF n (p!, . . .,p n , y/s,a 3 (ys)) S n (piy/s/s, . . . ,p n \fs/Tj . (4) 

Here there is a sum over the number n of final state partons. For a next-to-leading order 
calculation of a three-jet observable, n runs from 3 to 4. The cm. energy is y/S. There 
are integrations over dimensionless momenta pj for the final state partons. The true, di- 

mensionful momenta are Pj = pj \J S/ s where yfs = \Pj\- There is a delta function for 
momentum conservation. In place of the delta function for energy conservation, there is a 
function h(y/s) that satisfies / d^fs h(y/s) = 1. There is a conventional normalizing factor 
l/°o where a is the Born cross section for e + + e~ — > hadrons. 

The functions F n contain the perturbative matrix elements, calculated from the Feynman 
diagrams for the process. For instance, two such diagrams are depicted in Fig. 1. The first 
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FIG. 1: Two cuts of one of the Feynman diagrams that contribute to e + e — * hadrons. 

has four partons in the final state and contributes to F 4 , while the second has three partons 
in the final state and contributes to F3. In the case of diagrams with a virtual loop, as in 
the second diagram in Fig. 1, the function F n contains an integral over a dimensionless loop 
momentum I. It is an essential part of the numerical method used here that the integration 
over the loop momentum is performed numerically, at the same time as the integrations over 
the final state momenta are performed. 

At next-to-leading order, the functions F n contain either one or two powers of a s . We 
take the renormalization scale to be \i = c^^/S, where c M is a dimensionless parameter. In 
our numerical example at the end of this paper, we take c M = 1/6. 

The last ingredients in Eq. (4) are the measurement functions <S n , which are functions of 

the dimensionful final state momenta Pj\JS/ s. These functions are depicted in Fig. 1 by dots 
where parton lines cross the cut that represents the final state. The functions S n represent 
a three-jet observable in the sense that S2 = 0. They are invariant under interchanges of 
their arguments and are infrared safe in the sense that 

S n +l(Pl, XPn, (1 - X)Pn) = S n (Pl, ...,P n ) (5) 

for < A < 1, so that the measurement is unchanged if a parton splits into two partons 
moving in the same direction. 



III. ONE LOOP SELF-ENERGY DIAGRAMS 

In order to formulate parton showering properly, we need to understand how one parton 
splits into two at the one loop level in the Coulomb gauge. In this section, we analyze the 
perturbative calculation of a cut quark propagator that occurs as part of a larger graph. 
The analysis for gluons is essentially the same. 

To set the notation, consider first the contribution at order 

J[Born]= /^Trpo}, (6) 

where q is the three-momentum carried by the propagator, g M = (|<f |, q), and Rq denotes the 
factors associated with the rest of the graph and with the final state measurement function. 
In particular, Rq has two Dirac indices that match the Dirac indices on the $ and there is a 
trace over these indices. 

Next, consider the contribution from a cut self-energy insertion on the quark propagator 
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in question. The quark splits into a quark plus a gluon, giving a contribution of the form 



Here there is an integration over the loop momentum that describes the splitting. We use 
coordinates {q 2 , x, 0} for this loop momentum. We define these as follows. Let the gluon 
momentum be k + and the quark momentum after the splitting be k-, with k + + k- = q. 
We define q 2 and x by 



^/l + q 2 /\q\ 2 = (|* H 

2x-l = (\k + \-\k_\)/\q\. (8) 

Then < q 2 and < x < 1. The variables q 2 and x are constant on surfaces that are, respec- 
tively, ellipsoids and hyperboloids in loop momentum space. These surfaces are orthogonal 
to each other where they intersect. The virtuality of the quark-gluon pair is 

(\k + \ + \k-\) 2 -\q\ 2 = q 2 . (9) 

Finally, we let <fi be the azimuthal angle of k + in a coordinate frame in which the z axis lies 
along q and the x axis is defined arbitrarily. 

The function [a s /(2ir)]M g / q (q 2 , x, 0) in Eq. (7) is determined directly from the Feynman 
rules and is given as M q in Eq. (87) of Ref. [8]. 1 Its q 2 — > limit is simple: 

^M a / q {?,x,<l>)~i?±P g/q (x), (10) 

where P g / q (x) is the one loop Altarelli-Parisi kernel for q — » g, namely Cp[l + (1 — x) 2 ]/x. 
[We follow the notation of Ref. 4, in which P a /t,(x) denotes the Altarelli-Parisi kernel for 
b — > a but omitting regulation of the x — > 1 singularity, if there is such a singularity for that 
{a, b} combination.] 

The function R denotes the rest of the graph and the measurement function, as before. 
Because it includes the measurement function, it depends on x and as well as on q 2 . 
Because the measurement function is infrared safe, the q 2 — > limit of R is simple: 

i?(g 2 ,x,0) -> Rq as q 2 -> 0, (11) 

where R is the function appearing in the Born cross section, Eq. (6). 

Next, consider the contribution from virtual self-energy insertions on the quark propaga- 
tor in question. The quark splits into a quark plus a gluon that recombine, leaving a single 
quark line in the final state. The net contribution from the two virtual loop graphs is 

XMual] = / J| Tr { - f f /> jit I *«/.«".*> *A ■ < 12 > 



1 For the purposes of this paper, we thought it best to display the factor a s explicitly. 
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FIG. 2: Graphical representation of the three terms in Eq. (14). The third diagram is intended to 
represent the entire virtual contribution. 

Here, again, = (\q\,q) is the on-shell incoming momentum. The function V g / q is given 
in the Appendix. It equals the function V q calculated in Ref. [8] except for modifications to 
the terms that reflect renormalization. The most crucial feature of V g / q is its q 2 — > limit: 

Vg/g(q 2 ,X) -> Pg/ q (x) &S f ~> . (13) 

If we combine the three perturbative contributions, we have 



X[Born] + X[real] + J[virtual] = [ -% Tr(^ + [°° % C dx [ 

J 2 o I Jo q 2 Jo J- 



. . .. dx — 

2\q\ I Jo q 2 Jo J--k2tx 

^M g/q (q 2 } x,(f))R(q 2 ,x } (j)) - ^ V g/q (q 2 , x) ^R G 



.(14) 



These three terms are represented in Fig. 2. The real and virtual contributions are separately 
infrared divergent. However the integrands are added before the integration is performed 
numerically. Then the divergence cancels because of Eqs. (10), (13), and (11). (It is also 
important for the cancellation that the two terms match when q 2 — > and x — > at the 
same time, with q 2 /x fixed.) 



IV. ONE LEVEL SHOWER FOR A SINGLE PARTON 

So far, we see how a next-to-leading order numerical calculation works, but there is no 
parton shower. If we wanted to have the first splitting in a parton shower, we could follow 
the spirit of the parton shower Monte Carlo programs [1-3] and replace Eq. (14) by 

_ r . , f dq _, f r°° dq 2 f 1 , r w d(p a s . . , 2 lSr ,/-2 ,\ 
j [sh0 wer] = J m ^[l ^1 dxj_^- ^M„„(q ,xA)R(i ,*,*) 

*~°{-£Tl ld *%™™)\ (15) 

The quantity Xfshower] is represented graphically in Fig. 3. 

In Eq. (15), we keep the construction simple by maintaining a s evaluated at a fixed 
scale fi. We choose not to take the q 2 — > limit in any part of the function R(q 2 ,x,<f)) 
that describes the rest of the graph. The parton splits with the exact matrix element M. g / q , 
whereas a typical Monte Carlo program would use the small q 2 limit of M. g / q - As is standard, 
the exponential factor is interpreted as the probability for the parton not to split at virtuality 
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FIG. 3: Graphical representation of Xfshower] in Eq. 15. The pair of square vertices represents the 
matrix element Ai and the Sudakov exponential. 

greater than q 2 . In place of V g / q in the exponent, one might use the Altarelli-Parisi kernel 
supplemented by some theta functions that keep x away from and 1, but we use the exact 
function V g / q that occurs in the virtual corrections. 

We see here that parton shower calculations as exemplified in Eq. (15) are very different 
from fixed order calculations as exemplified in Eq. (14). In the fixed order calculation, 
the virtual graph acts as a subtraction to the real splitting contribution. In the shower 
calculation, the virtual graph becomes a multiplicative modification to the splitting matrix 
element. 

The perturbative representation in Eq. (14) is the same as the shower representation in 
Eq. (15) up to corrections of order a 2 : 



J[Born] + X[real] + Zfvirtual] = X[shower] x (l + 0(a\ 



(16) 



The general idea incorporated in this subtraction to multiplication theorem is a standard part 
of the strategy behind shower Monte Carlo programs, 2 but we have not found the theorem 
stated in just this way, so we prove it here. First, we write Xfshower] as 



with 



T\ [shower] 



dq 

W\ 



Tr 



Xfshower] = X\ [shower] + X 2 [shower] 
dq 2 n 



(17) 



o q 2 



dx 



^2tt 



exp - 



oo dl 2 
* I 2 



2tt 



dz ^-Vg/ q (P, Z) 



x 



a 



^M g/q (q ,x,(/>)R(q ,x, 



Re 



(18) 



and 

1 2 [shower] 



dq 





poo 


- Trj 












4 


dx — 




Jo 2tt 



dq 2 



» dl 2 
J 2 JU 



a-. 



dZ ^K Valq ^ 



(19) 



Let us find the lowest order (up to and including order ct s x R) terms in the expansions 
of these integrals. The integral 1\ [shower] is simple. We take the order ot° s term in the 



See, for example, Eq. (B.15) in Appendix B of Ref. [11]. 
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exponential to obtain 



X l[s hower] = j^\[ d jrj^f_2 

' +0(a 2 s x R). (20) 



cx 

^M g/g {q 2 ,x,<i))R{q 2 ,x,<i)) - V g/q {q 2 , x) iR 

The integral Z 2 [shower] is more subtle. We cannot simply expand under the integral since 
the integrals of each term would be divergent at small q 2 . Instead, we note that we have 
the integral over q 2 of a total derivative, so that the integral is just the difference of the 
integrand between the integration end points: 

J 2 [shower] = Tr {4 R } [F(oo) - F(0)], (21) 



where 



Now the function V g / q (J 2 ,z), which is given in the Appendix, has the important property 
that it is proportional to 1/7 2 for I 2 — > oo. This insures that the integral in the exponent 
is convergent for I 2 — > oo, so that F(q 2 ) — > exp(— 0) = 1 as q 2 — > oo. On the other hand, 
Vg/ q (J 2 ,z) approaches P g / q (z) as I 2 — > 0. This makes the P-integral diverge as q 2 — > 0. In 
fact, one gets two powers of log(g 2 ) for q 2 — > 0, one directly from the P-integral and one 
from a divergence at z — > in the z-integral that develops as V g / q (P, z) gets closer and closer 
to Pg/ q (z). The net effect is that the exponent behaves like a positive constant times \og 2 (q 2 ) 
for q 2 0. Thus F(q 2 ) — > exp(— oo) = as q 2 — > 0. The result for Z 2 [shower] is then 

J 2 [shower] = f^-TrURo}. (23) 

Adding the results (20) and (23), we have the result in Eq. (14) plus 0(a 2 x R). This proves 
the theorem. 

There is an analogous result for the gluon propagator. The trace is over Lorentz indices 
instead of Dirac indices. The functions M. g / q and V g / q are replaced by functions M- g / g + 
N F A4 q / g given 3 in Ref. [8] and V g / g + N F V q / g given in the Appendix. 

It should be evident that one can use the subtraction to multiplication theorem to add 
parton showers to the next-to-leading order calculation of three-jet cross sections, while 
keeping unchanged the first two terms in the perturbative expansion of the result. Indeed, 
there is more than one way to do this. Which way is best can and should be debated, but 
in this paper we confine ourselves to mapping out one way. 

We will refer to the first splittings of the partons emerging from a Born graph as primary 
splittings. The basic idea for the primary splittings is this. Start with a computer program 
that does calculations at next-to-leading order in the Coulomb gauge. In each cut a] graph, 
replace X[Born] by Xfshower] for each cut propagator. Then take the cut order a 2 graphs 
that have a cut self-energy subgraph plus its corresponding virtual correction. In each such 



3 To be precise, [a s /(2Tv)](M g / g + N F M q / g ) is given in Eq. (37) of Ref. [8] as M g . 
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graph, replace X[real] + Xfvirtual] by 0. (There are some complications that arise when 
applying this basic idea. We address the complications in subsequent sections.) 

Once one understands this basic idea, one can envisage a lot of refinements. Here we will 
pursue one possible refinement. In X[shower] we choose a number Ay of order 1 and divide 
the integration region into a low virtuality part, < q 2 < Ay<f 2 , and a high virtuality part 
Ay<f 2 < q 2 < oo. In the high virtuality term, we expand the Sudakov exponential in powers 
of a s and keep only the order term. This gives 

Xfshower] = (x[shower, Ay] + J[real, Ay]) x (l + C(« 2 )) (24) 

where 

_ r . , , f dq _ f r x v<l 2 dq 2 r l , / >7r d(p a s . . , , , N „/_ 2 

Z[shower, V] = / ^ / o * / ' *' W« • *) 

and 

X[real,A y ] = / ^ Tr { f° *f f ^ ^-M s / 9 (g 2 , x, 0) i?(g 2 , x, 0)) . (26) 

J 2|g| {J\ v q 2 q Jo J-ttZtt 2tt ' J 

Then the refinement of the basic idea for the primary splittings is as follows. Start with 
a computer program that does calculations at next-to-leading order in the Coulomb gauge. 
In each cut a] graph, replace X[Born] by Xfshower, Ay] for each cut propagator. Then take 
the cut order a 2 graphs that have a cut self-energy subgraph plus its corresponding virtual 
correction. In each such graph, replace X[real] + X[virtual] by X[real, Ay] . This approach is 
a little more complicated than the simple approach we explained first, which amounts to 
choosing Ay = oo. 

Why would one choose Ay < oo? Unless one has a NNLO calculation to compare to, one 
cannot know for sure which approach is best, since the approaches differ only at order af +2 . 
However, the Ay < oo approach has the advantage that we apply the formula with Sudakov 
suppression of small virtuality splitting in the region of small or moderate virtuality, while 
we use ordinary fixed order perturbation theory in the large virtuality region where "parton 
shower" physics is not relevant. 

Another possibility would be to modify X[shower] by replacing Ai by a new function 
M.' and by replacing V by a new function V. The new functions could be anything we 
like provided that they have the same small q 2 — > limits as the original functions. They 
could, for example, be the splitting functions in a "standard" showering scheme based on 
virtuality as the evolution variable. Then, in the order a 2 graphs, Ai would be replaced 
by M. — Ai' and V would be replaced by V — V. Thus the shower functions Ai' and V' 
would act as collinear subtraction terms in the order a 2 graphs. This is similar to what is 
done for initial state collinear singularities in [11, 12]. The treatment discussed above with a 
virtuality cutoff Ay < oo is a special case of this more general possibility. We do not pursue 
the general Ai', V possibility here. 

V. PRIMARY SPLITTINGS FOR A WHOLE GRAPH 

Let us see how this algorithm works in a particular example, taking Ay = oo for simplicity. 
Begin with a cut Born (i.e. order a]) graph with the topology shown on the left in Fig. 4. 
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FIG. 4: A cut Born graph. We depict here only the topology of the graph, leaving unspecified 
which lines are quarks and which are gluons. On the right, we show the graph after replacing the 
cut propagators by X[shower], which is represented by a circle containing an S. 




FIG. 5: Six order a 2 s cut Feynman graphs. (For a graph with a self-energy insertion with one 
adjacent propagator cut, we understand the complete virtual correction.) These graphs are the 
order o? s term in the expansion of the graph in Fig. 4. 

The primary splitting algorithm is to replace X[Born] by X[shower] for each of the three 
parton lines crossing the cut. This is illustrated in the right hand side of Fig. 4. 

When this expression is expanded in powers of a s , it gives the Born graph back along 
with the six order a 2 graphs depicted in Fig. 5. There are further contributions that are of 
order oP s and higher. Since we get the graphs in Fig. 5 from the expansion of the right hand 
diagram in Fig. 4, we delete these contributions from the order a 2 graphs. 

VI. PRIMARY SPLITTINGS FOR ANOTHER GRAPH 





FIG. 6: A cut Born graph and its replacement according to the level 1 showering algorithm. 
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There is one more topology for a cut Born graph that contributes to three-jet cross 
sections. This topology is exemplified by the graph shown as the left hand diagram in 
Fig. 6. Again, we keep this analysis as simple as possible by treating the case that the 
virtuality cut is Ay = oo. 

To treat this graph, we make use of a technique that helps keep two-jet physics from 
affecting the three-jet cross section. Let k\ be the momentum carried by the top propagator 
in Fig. 6 and let k 2 be the momentum carried by the bottom propagator before it splits into 
two propagators. We multiply the graph by 

l = F 1 (k 1 ,k 2 ) + F 2 (k 1 ,k 2 ), (27) 

where 

*> y = JW' (28) 

We separate the two terms and treat them separately, considering the factor Fj as part of 
the graph. We can indicate which term is which by drawing a line through the propagator 
that carries momentum kj when we have multiplied by Fj. The line indicates that the 
corresponding propagator is blocked from going on-shell because of the factor {kj) 2 . For 
the graph under consideration, only the F 2 term contributes, as indicated by the equality 
between the left and middle graphs of Fig 6. 

Now we adopt as the splitting algorithm for this case the replacement of each cut prop- 
agator by X[shower], as indicated by the right hand graph in Fig 6. It is significant that we 
do this only after inserting the factor F 2 . 

The role of the factor Fj is easy to understand on an intuitive basis. It is possible for 
the graph on the left in Fig. 6 to generate a two-jet event in which the gluon is nearly 
parallel to the antiquark at the bottom of the graph. When each parton is allowed to split, 
one could get a high virtuality splitting from the quark at the top of the graph, creating a 
three-jet event. Of course, we want to generate three-jet events, but our calculation is based 
on generating three-jet events directly from the Born graphs. The function Fj enforces the 
requirement that the high virtuality splitting be the one in the Born graph. 

Consider now the eight order a 2 cut Feynman graphs illustrated in Fig. 7, where we have 
inserted factors of 1 = Fi + F 2 . For all except the first graph shown, only the F 2 term 
survives, as indicated. For the first graph, there is a corresponding graph with a factor F±, 
which is not shown and would be grouped with a different set of graphs. 

Now notice that when the right hand graph of Fig. 6 is expanded in powers of a s , it 
reproduces the Born graph with which we started plus the graphs of Fig. 7 plus corrections 
that are of order or higher. Thus we take the splitting algorithm for this case to be that 
the graphs of Fig. 7 are deleted from the calculation. 



VII. THE NEXT-TO-LEADING ORDER GRAPHS 

We have seen how to insert a one level parton splitting into the Born graphs. We have 
also seen that certain of the order a 2 cut graphs are thereby generated, so that these graphs 
are to be deleted from the calculation. In the case that we have chosen \y < oo, then the 
graphs that would have been deleted are retained, but with a virtuality cut q 2 > Xyq 2 on 
the appropriate cut self-energy subdiagram. Other order a 2 cut graphs are not generated 
from the Born graphs and need to be added to the Born contributions. 
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FIG. 7: Eight order a 2 cut Feynman graphs. 



VIII. SECONDARY SPLITTINGS 



We have replaced each cut propagator in the Born diagrams with one step of a parton 
shower and have called the parton splittings generated in this way primary splittings. Now 
we would like to let these partons split further with secondary splittings. For each of the 
parton lines entering the final state, we replace X for that line by a simplified version of the 
splitting probability. For example, if the parton is a quark, we use 



X[shower-0] = f TrJ J ^- [ dx [ ^ ^ S g/q (q 2 , x) $ R(q 2 , x, 
J 2\q\ | Jo q z Jo J-irzn Z7r 



x exp — 



dq 
Q 

°° dl 2 
2 "F 



a 2 rl 



n d(p a s 



2tt 



dz 7^-S g / q {P,z) 



(29) 



where the Sudakov exponent S g / q is a suitable simplified version of V g / q . Here the Dirac 
structure, ^, is kept exactly the same as it was in X[Born] for this line. We thus take q° = \q\ 
in ^. Furthermore, we take q° = \q\ in the rest of the graph. When this expression is 
expanded in powers of a s , it gives the original graph back along with some order and yet 
higher order corrections. 

In the construction so far, we have used a coupling a s evaluated at a fixed scale, /i. As 
a refinement, we follow standard practice by using for the showering in Eq. (29) a coupling 
evaluated at a running scale. We take the scale of a s in the main line of the equation to 
be [x(l — x)q 2 S/ sq] 1 I 2 and in the exponent to be [z(l — z)PS/so\ 1 ^ 2 . Recall that we use a 
numerical integration program that employs dimensionless momenta as integration variables. 
The factor x(l — x)q 2 is the dimensionless squared transverse momentum that is associated 
with the splitting. A factor S/s would convert to the physical momenta. We use instead 
S/sq, where y/s^ — is the sum of the absolute values of the dimensionless momenta in 

the perturbative diagram before splittings. Of course, when we introduce a running coupling 
inside an integral as in Eq. (15), we create an ambiguity about how to evaluate a s (k) when 
k < Aq CD . We resolve the ambiguity by freezing a s (k), so that it ceases to grow with 
decreasing k as soon as it reaches the value a s — 1. 

The function S g / q (q 2 , x) should be fairly simple, its integral in the Sudakov exponent 
should be finite, and it should approach the Altarelli-Parisi function P g / q (x) as q 2 — > 0. The 
reader may choose his or her favorite function. We have tried the following choice, which 
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satisfies these criteria, for the purposes of this paper, 



S g/q (q 2 ,x) = P g/q (x) 9(q 2 < x(l - x)k 2 ). (30) 

Here k 2 serves as a cutoff parameter for the splitting. The theta function may be regarded 
as a simplified version of the factors 4x(l — x)/[q 2 /q 2 +4x(l —x)] that occur in the functions 
V that represent the virtual self-energy diagrams and appear in the Sudakov exponent of 
the primary splittings. (See the Appendix for the functions V '.) We examine what to choose 
for k below. 

We can now move on to generate parton splittings within the previous parton splittings. 
We do this indefinitely except that, when q 2 for a splitting is smaller than some prede- 
termined cutoff, we regard the splitting as having gone too far. We then replace the two 
daughters by the mother, canceling the splitting. The mother parton remains in the final 
state and no further splitting of it is attempted. 

Since the mother parton in the current splitting was created in a previous splitting, one 
can choose k 2 so as to insure angular ordering between the two splittings. For a nearly 
collinear splitting with parameters {q 2 ,x, \ q\} of a mother parton with momentum q, the 
angle 9 between the two daughter partons is given approximately by 

02 = n ~ v* - (31) 

x(l — x)q z 

Suppose that the previous splitting had parameters {p 2 , y, \q\/y}, with y being the momen- 
tum fraction of the daughter parton that becomes the mother parton for the second splitting. 
Then the angle between the two partons in the first splitting can be estimated as 



9q min 



yp 2 _ ' 



(i -vW 



(32) 



The first expression here uses the small angle approximation, while the second provides a 
rough upper bound for the case that the angle is not small. Thus if we take 



k 2 = min 



VP 2 ^ 2 2 

" j Q j ^mother 



y) 



(33) 



the probability for the second splitting is zero when the second angle is greater than the 
first. Here ^mother is the k parameter that was used for generating the splitting that gave 
the mother parton in the current splitting. 

We have not yet specified what the splitting scale k 2 should be for the first secondary 
splitting, the splitting of one of the partons that was generated in a primary splitting. We 
should be a bit careful because the function R used in the proof in Sees. Ill and IV now 
includes the measurement function convoluted with the functions describing the secondary 
splittings. In the notation of the preceding paragraph, R(p 2 , y, (f>) should approach i? when 
p 2 — > 0. This will happen if k: 2 — ^ (so that the secondary jets become narrow) when 
p 2 — > 0. For this reason, we take 



k 2 = min 



yp 2 ^2 -2 

Q ,c K p 



.(i-yY 

where c K is a constant, taken to be 4 in the numerical example studied in Sec. X. 



(34) 
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IX. SINGULARITIES OF THE NEXT-TO-LEADING ORDER GRAPHS 



In the next-to-leading order calculation without showers, the order a 2 cut graphs contain 
divergences arising from soft and collinear parton configurations. These divergences cancel 
between different cuts of a given graph, that is between real and virtual parton contributions. 

In the Coulomb gauge, the collinear divergences arise from self-energy graphs. When 
we incorporate parton showers the collinear divergences disappear. If we use Eq. (16), we 
simply eliminate the self-energy graphs associated with final state particles, and with them 
we eliminate the collinear divergences. With the small modification in Eq. (24), we retain 
part of the real splitting self-energy graph in which the virtuality q 2 is bigger than XyQ 2 ■ 
This still eliminates the collinear divergence. As mentioned in Sec. IV, we could adopt 
suitably behaved alternative functions M.' and V to replace M. and V in the splitting from 
a Born level graph. Then the order a 2 s real and virtual self-energy graphs would still be 
present, but Ai would be replaced by Ai — M! in the real graphs and V would be replaced 
by V — V in the virtual graphs. Thus the shower functions Ai' and V would act as collinear 
subtraction terms in the real and virtual order a 2 graphs. In this case also, the collinear 
divergences are eliminated separately from the real and virtual graphs. 

This leaves only the soft divergences in the order a 2 graphs. The soft gluon singularities 
are an essential part of QCD. Their structure is quite different from that of the collinear 
singularities. Their treatment will be the subject of another paper [9]. Once the soft singu- 
larities are accounted for, one can let each parton from an order a 2 graph create a shower. 
Meanwhile, the treatment given here, without showers from the order a 2 graphs, yields a 
divergence-free calculation because the soft gluon divergences cancel between different cuts 
of a given graph. We will exhibit some results from this calculation in the following section. 



X. A NUMERICAL TEST 

In this section, we provide a numerical test of the algorithm presented here. We have 
computed the thrust distribution da/dt for thrust t = 0.86, in the middle of the three-jet 
region. We compute the thrust distribution using the algorithm described in this paper. 
That is, for self-energy graphs in the limited virtuality region, we replace " Born x [1 + real — 
\virtual\]" by "Born x real exp(— \virtual\) Secondary showering from the Born graphs is 
also included. On the other hand, in this test there are no showers from the order ctf +1 
graphs. 

The idea is to compare the thrust distribution thus calculated to the same distribution 
calculated with a straightforward NLO computation with no showers. According to the 
analysis presented above, the difference between the two calculations should be of order 
af +2 . We divide the difference by the NLO result, forming 

(NLO-shower) - NLO 
R ~ NLO • (35) 

The ratio R should have a perturbative expansion that begins at order a 2 . We can test 
this by plotting the ratio against a 2 . We expect to see a curve that approximates a straight 
line through zero for small a 2 . For comparison, we exhibit also the ratio with the af +l 
corrections omitted, 

(LO-shower) - NLO 
^ LO = NLO • (36) 
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FIG. 8: Comparison of the NLO calculation with showers partially added as described in this 
paper to a pure NLO calculation using [10]. We plot the ratio R defined in Eq. (35) for the 
thrust distribution at thrust equal 0.86. Also shown is the ratio i?LO; defined in Eq. (36), in 
which the order af +1 correction terms are omitted from the calculation. The cm. energy is 
y/~S = Mz and the renormalization scale is chosen to be fi = \/S/6. These ratios are calculated 
for a s (M z ) 2 = {0.25,1,2,3,4} x (0.118) 2 and plotted versus a s (M z ) 2 /(0.118) 2 . 

The ratio i?Lo should have a perturbative expansion that begins at order af Thus we expect 
to see a curve proportional to the square root function [a 2 ] 1 ^ 2 for small a 2 . (The size of the 
coefficient of [a 2 ] 1 / 2 in R LO does not have any great significance, since it is quite sensitive 
to the choice of renormalization scale.) A comparison of the two curves, R and R^o, shows 
whether the af +1 corrections, which are lacking in R^o, are doing their job. 

The results of this test are shown in Fig. 8 for v^S* = Mz and for a range of a s [Mz) 
from (1/2) x 0.118 to 2 x 0.118, where 0.118 represents something close to the physical 
value for a s (M z ). We choose the renormalization scale to be \i = y/S/6. One should note 
that 2 x 0.118 amounts to quite a large a s in this calculation since a s [Mz) = 0.236 gives 
a 8 (fj) = 0.586. 

We see the expected shape of the i?LO curve. We also see that the R curve lies closer 
to zero, indicating that the af +1 corrections are acting in the right direction. The more 
exacting test is to examine whether the R curve approaches a straight line through the origin 
as a\ — > 0. Within the errors, it does. However, the slope of the straight line is quite small. 
That is presumably an accident of the choice of parameters in the program. It remains for 
the future to examine the effects of various parameter choices. 

In Fig. 9, we show the same comparison, but this time for t = 0.71. This is near the 
value t = 2/3 that marks the far end of the three-jet region, with the three partons in what 
is sometimes called the Mercedes configuration. In this case, the R curve again approaches 
a straight line through the origin as a 2 — > 0, again with quite a small slope. 

In Fig. 10, we show the same comparison, but this time for t = 0.95. This is near the 
two-jet limit at t — 1. For a s = 0.118, one would normally not use a calculation that did 
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FIG. 9: Comparison of the NLO calculation with showers partially added to a pure NLO calculation 
for t = 0.71. The notation is as in Fig. 8. 

not include a summation of logs of 1 — t for t this close to 1, since log(0.05) 2 ~ 9. Thus we 
would not recommend using the code discussed in this paper for a comparison to data this 
near to the two-jet limit. Nevertheless, we can still test for the absence of an a] term in R. 
Looking at the graph, we see that, within the errors, there is no evidence for a nonzero a I 
term in R. We expect an a 2 term, but the coefficient of a 2 appears to be quite small. On 
the other hand, it appears that some of the higher order terms are quite substantial. 

XI. CONCLUSIONS 

We have presented a method for adding parton showers to a next-to-leading order calcu- 
lation in QCD. Specifically, we treat the final state collinear singularities, using the example 
of e + + ^3 jets. Our intention is that a program like ours would be used in conjunction 
with a standard Monte Carlo event generator [1-3]. The complete calculation would then 
act as an event generator (with positive and negative weights for the events) in which the 
final states consist of hadrons generated from parton showers. Recall that we defined pri- 
mary parton splittings to be those of the partons emerging from an order af graph. There 
are also secondary splittings: the further splittings of the daughters from the primary split- 
tings as well as the splittings of partons emerging from an order af +1 graph and of their 
daughters. It is the primary splittings that need to be matched to the NLO calculation. 
The secondary splittings affect an infrared safe observable only at order af +2 . Thus the 
partonic events from our program after the primary splittings could be handed to a Monte 
Carlo event generator, which would perform the secondary splittings and the hadronization. 
Alternatively, our program can perform the secondary splittings, leaving hadronization for 
the Monte Carlo event generators. For the purpose of implementing such a scheme, it is 
significant that the weight functions do not have singularities when two partons are collinear 
or one is soft. This is in contrast to pure NLO calculations, which depend on cancellations 
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FIG. 10: Comparison of the NLO calculation with showers partially added to a pure NLO calcu- 
lation for t = 0.95. The notation is as in Fig. 8. 



between different events when an infrared safe observable is calculated. 

The method presented here is incomplete in that we leave for a companion publication 
[9] the treatment of the soft gluon singularities of QCD. The code described in this paper is 
available at [10]. A treatment of soft gluon effects, to be described in [9], is included in the 
code. Parton showers are included, but the requisite data structures and color flow informa- 
tion that would be necessary to provide an interface with the showering and hadronization 
models of standard event generator Monte Carlo programs are not yet available. 
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APPENDIX: VIRTUAL SPLITTING FUNCTIONS 



In this appendix, we examine the functions V a /b(q 2 ,x) that specify the virtual parton 
self-energy graphs in the Coulomb gauge. 

For quark splitting, the function that we use is 

^(f^) = {12x(l-x)-l} 
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-r- { 20a; ( 1 - *) - 3 + log^V/^Ld) [12^(1 -x)-i]} 

AC 

+-^x(l-x){5- Ux(l-x)} 

x(l-x){l- 6x(l -x) + S[x(l - x)] 2 } 

^x 

8C p , „ 

-Djx} 2X - 1)1(1 - X) 
1 fiC 

~^m {2x ~ l)x{l ~ x){l ~ 2x{l ~ x)h (A - 1} 

where /x is the MS renormalization scale, 

A*Ld = min(/i 2 , |g | 2 ), (A.2) 

and 

^ = q 2 /\q\ 2 + ±x(l-x), 

A = ?Vk1 2 + l- (A.3) 

The function V g / q (q 2 ,x) has been modified from the function V q given in Ref. [8]. We have 
replaced the two terms 

2 e- 3 g 2 //i 2 + l F e- 5 / 3 g 2 //i 2 + l 1 
by the first two terms in Eq. (A.l). The contributions of these terms to the integral 

-wV g f q {f,x) (A.5) 



r°° a. 
U 1 



are the same in the limit of small q 2 up to terms that vanish as q 2 — > (and are non- 
singular as functions of x). Thus the two expressions are equivalent when inserted into the 
perturbative formula (14). 

The terms that we have modified contain the effects of MS renormalization. In the 
original form, they cut off the integral (A.5) at I 2 of order /x 2 . That is a sensible way to 
express perturbative renormalization. However, in our present application the integral 

roc tfp- rl 

j- 2 jr i d^ g/q (P,x) (A.6) 

appears in the Sudakov exponent for shower generation. In the case that \q\ 2 <C /x 2 , the 
former terms in Eq. (A. 4) would have caused the Sudakov exponential to deviate from 1 for 
q 2 much greater than \q | 2 . The effect of this deviation on the calculated cross section is of 
order a 2 x R and thus is, in principle, negligible. However this effect appears to us to be 
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quite unphysical. The revised form of V g / q is preferable in that it confines the contribution 
of the renormalization dependent terms to q 2 of order \q | 2 or smaller. 

One may also note that the integral (A. 6) is positive for all q 2 and increases as q 2 decreases. 
This was not the case with the form of V g / q before the replacement. 

For gluon splitting to a quark and antiquark, we use 

V q/9 {q\x) = {1 " 2 < l - *)> ( A - 7 ) 



^2 3-{ 2 -!^^) + log(/i 2 //iLd)[l - 2x(l - x)}} . 



For gluon splitting to two gluons, we use 

V 9/9 {q\x) = 5* {3x(l - x) - 1} 



+t4 75- {7*(1 " ^ " 2 + log^/^Ld)^!! -*)-!]} 

it rmod 



+^^x(l-x){l-2x(l-x)} 

+ x(l-x){l- Ax(l -x)+ A[x(l - x)] 2 } . (A.8) 

x 

Here the definitions (A. 2) and (A. 3) apply again. For V q / g , we have replaced 

1 1 X(1 ~ X) (A.9) 



2 e~ 2 q 2 /fi 2 + 1 e~^q 2 /n 2 + 1 
given in Ref. [8] by the terms in Eq. (A. 7). For V g / g , we have replaced 

r 1 i r x{i-x) x(l-x) 

UA e- 2 q 2 /^ 2 + 1 + ^ A e^q 2 /^ 2 + 1 + ^ V^g 2 /^ 2 + 1 1 UJ 

given in Ref. [8] by the first two terms in Eq. (A.8). The reasoning behind this modification 
is the same as in the case of quark splitting. 

The Sudakov exponent for gluon splitting is proportional to 



r°° dl 2 r 1 

J 5* I 2 Jo 



dx 



V g/g (l 2 ,x) + N F V q/9 (l 2 ,x) . (A.ll) 



This quantity is positive for all q 2 and increases as q 2 decreases provided that N F > 3C A /2. 
Thus this form will work fine for the main cases of phenomenological interest, Ca = 3 with 
5 or 6 quark flavors. Further modifications would be necessary for other cases. We have not 
pursued this topic further because we expect that the V functions will be different anyway 
in future versions of the algorithms presented in this paper. 
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